## This code creates Table 3, Figure 8 (a)-(d).

neighbor_dist_cat <- seq(10,100,10)

# Table 3, 1912 results

data <- colonial_towns_shp[colonial_towns_shp$neighbor_dist <= 50 & colonial_towns_shp$year==1912,]

force_specs <- c("total_force_OLS_BAGO_NoGEO","total_force_OLS_BAGO_NoFE","total_force_OLS_BAGO_FULL","m_force_OLS_BAGO_NoGEO","m_force_OLS_BAGO_NoFE","m_force_OLS_BAGO_FULL")

total_force_OLS_BAGO_NoGEO <- conleyreg(total_force_ih~replaced+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
total_force_OLS_BAGO_NoFE <- conleyreg(total_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
total_force_OLS_BAGO_FULL <- conleyreg(total_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)

m_force_OLS_BAGO_NoGEO <- conleyreg(m_force_ih~replaced+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_force_OLS_BAGO_NoFE <- conleyreg(m_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_force_OLS_BAGO_FULL <- conleyreg(m_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)

force_table <- matrix(nrow = 7,ncol=7)

for(i in 1:length(force_specs)){
  force_table[1,i+1] <- force_specs[i]
  model <- get(force_specs[i])
  coef <- round(model$coefficients[2],digits=4)
  se <- round(model$coefficients[2,2],digits=4)
  force_table[2,i+1] <- coef
  force_table[3,i+1] <- paste0("(",se,")")
  force_table[4,i+1] <- paste0("[",round(coef-1.96*se,digits=4),",")
  force_table[5,i+1] <- paste0(round(coef+1.96*se,digits=4),"]")
  force_table[6,i+1] <- nobs(model)
  force_table[7,i+1] <- round(model$adj.r.squared,digits=4)
}

xtable(force_table)

# Figure 8 (a)

total_force_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(total_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1912 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  total_force_BAGO[j,"coefficient"] <- model[2,1]
  total_force_BAGO[j,"se"] <- model[2,2]
  total_force_BAGO[j,"upper"] <- total_force_BAGO[j,"coefficient"] + (1.96*total_force_BAGO[j,"se"])
  total_force_BAGO[j,"lower"] <- total_force_BAGO[j,"coefficient"] - (1.96*total_force_BAGO[j,"se"])
  total_force_BAGO[j,"upper_10pc"] <- total_force_BAGO[j,"coefficient"] + (1.645*total_force_BAGO[j,"se"])
  total_force_BAGO[j,"lower_10pc"] <- total_force_BAGO[j,"coefficient"] - (1.645*total_force_BAGO[j,"se"])
  total_force_BAGO[j,"pch"] <- ifelse(sign(total_force_BAGO[j,"upper"])==sign(total_force_BAGO[j,"lower"]),"black",ifelse(sign(total_force_BAGO[j,"upper_10pc"])==sign(total_force_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/total_force_BAGO_1912.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-2,-1.6,-1.2,-0.8,-0.4,0,0.4,0.8,1.2,1.6,2))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,total_force_BAGO$lower/4,x.axis,total_force_BAGO$upper/4)
points(x.axis,total_force_BAGO$coefficient/4, pch=21, bg=total_force_BAGO$pch)
abline(h=0, lty=3)
dev.off()

# Figure 8 (b)

m_force_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(m_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1912 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  m_force_BAGO[j,"coefficient"] <- model[2,1]
  m_force_BAGO[j,"se"] <- model[2,2]
  m_force_BAGO[j,"upper"] <- m_force_BAGO[j,"coefficient"] + (1.96*m_force_BAGO[j,"se"])
  m_force_BAGO[j,"lower"] <- m_force_BAGO[j,"coefficient"] - (1.96*m_force_BAGO[j,"se"])
  m_force_BAGO[j,"upper_10pc"] <- m_force_BAGO[j,"coefficient"] + (1.645*m_force_BAGO[j,"se"])
  m_force_BAGO[j,"lower_10pc"] <- m_force_BAGO[j,"coefficient"] - (1.645*m_force_BAGO[j,"se"])
  m_force_BAGO[j,"pch"] <- ifelse(sign(m_force_BAGO[j,"upper"])==sign(m_force_BAGO[j,"lower"]),"black",ifelse(sign(m_force_BAGO[j,"upper_10pc"])==sign(m_force_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/m_force_BAGO_1912.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-2,-1.6,-1.2,-0.8,-0.4,0,0.4,0.8,1.2,1.6,2))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,m_force_BAGO$lower/4,x.axis,m_force_BAGO$upper/4)
points(x.axis,m_force_BAGO$coefficient/4, pch=21, bg=m_force_BAGO$pch)
abline(h=0, lty=3)
dev.off()

# Table 3, 1924 results

data <- colonial_towns_shp[colonial_towns_shp$neighbor_dist <= 50 & colonial_towns_shp$year==1924,]

force_specs <- c("total_force_OLS_BAGO_NoGEO","total_force_OLS_BAGO_NoFE","total_force_OLS_BAGO_FULL","m_force_OLS_BAGO_NoGEO","m_force_OLS_BAGO_NoFE","m_force_OLS_BAGO_FULL")

total_force_OLS_BAGO_NoGEO <- conleyreg(total_force_ih~replaced+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
total_force_OLS_BAGO_NoFE <- conleyreg(total_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
total_force_OLS_BAGO_FULL <- conleyreg(total_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)

m_force_OLS_BAGO_NoGEO <- conleyreg(m_force_ih~replaced+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_force_OLS_BAGO_NoFE <- conleyreg(m_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_force_OLS_BAGO_FULL <- conleyreg(m_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)

force_table <- matrix(nrow = 7,ncol=7)

for(i in 1:length(force_specs)){
  force_table[1,i+1] <- force_specs[i]
  model <- get(force_specs[i])
  coef <- round(model$coefficients[2],digits=4)
  se <- round(model$coefficients[2,2],digits=4)
  force_table[2,i+1] <- coef
  force_table[3,i+1] <- paste0("(",se,")")
  force_table[4,i+1] <- paste0("[",round(coef-1.96*se,digits=4),",")
  force_table[5,i+1] <- paste0(round(coef+1.96*se,digits=4),"]")
  force_table[6,i+1] <- nobs(model)
  force_table[7,i+1] <- round(model$adj.r.squared,digits=4)
}

xtable(force_table)

# Figure 8 (c)

total_force_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(total_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1924 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  total_force_BAGO[j,"coefficient"] <- model[2,1]
  total_force_BAGO[j,"se"] <- model[2,2]
  total_force_BAGO[j,"upper"] <- total_force_BAGO[j,"coefficient"] + (1.96*total_force_BAGO[j,"se"])
  total_force_BAGO[j,"lower"] <- total_force_BAGO[j,"coefficient"] - (1.96*total_force_BAGO[j,"se"])
  total_force_BAGO[j,"upper_10pc"] <- total_force_BAGO[j,"coefficient"] + (1.645*total_force_BAGO[j,"se"])
  total_force_BAGO[j,"lower_10pc"] <- total_force_BAGO[j,"coefficient"] - (1.645*total_force_BAGO[j,"se"])
  total_force_BAGO[j,"pch"] <- ifelse(sign(total_force_BAGO[j,"upper"])==sign(total_force_BAGO[j,"lower"]),"black",ifelse(sign(total_force_BAGO[j,"upper_10pc"])==sign(total_force_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/total_force_BAGO_1924.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-2,-1.6,-1.2,-0.8,-0.4,0,0.4,0.8,1.2,1.6,2))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,total_force_BAGO$lower/4,x.axis,total_force_BAGO$upper/4)
points(x.axis,total_force_BAGO$coefficient/4, pch=21, bg=total_force_BAGO$pch)
abline(h=0, lty=3)
dev.off()

# Figure 8 (d)

m_force_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(m_force_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1924 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  m_force_BAGO[j,"coefficient"] <- model[2,1]
  m_force_BAGO[j,"se"] <- model[2,2]
  m_force_BAGO[j,"upper"] <- m_force_BAGO[j,"coefficient"] + (1.96*m_force_BAGO[j,"se"])
  m_force_BAGO[j,"lower"] <- m_force_BAGO[j,"coefficient"] - (1.96*m_force_BAGO[j,"se"])
  m_force_BAGO[j,"upper_10pc"] <- m_force_BAGO[j,"coefficient"] + (1.645*m_force_BAGO[j,"se"])
  m_force_BAGO[j,"lower_10pc"] <- m_force_BAGO[j,"coefficient"] - (1.645*m_force_BAGO[j,"se"])
  m_force_BAGO[j,"pch"] <- ifelse(sign(m_force_BAGO[j,"upper"])==sign(m_force_BAGO[j,"lower"]),"black",ifelse(sign(m_force_BAGO[j,"upper_10pc"])==sign(m_force_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/m_force_BAGO_1924.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-2,-1.6,-1.2,-0.8,-0.4,0,0.4,0.8,1.2,1.6,2))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,m_force_BAGO$lower/4,x.axis,m_force_BAGO$upper/4)
points(x.axis,m_force_BAGO$coefficient/4, pch=21, bg=m_force_BAGO$pch)
abline(h=0, lty=3)
dev.off()

